Problema 1¶
Se tiene un oscilador harmónico en dos dimensiones, $(x,y)$ en coordenadas cartesianas,en el que la fuerza está dada por $\vec{F}=-k\vec{r}$ y se tiene una fricción dada por $\vec{f}_{fric}=\begin{cases} -\gamma\vec{v} &, \mbox{si } |\vec{v}|<1 \\ -\mu\ v^{3/2} \hat{v} &, \mbox{si } |\vec{v}|>1 \end{cases}$. Si $m=1, \gamma=0.1 \mbox{ y } \mu=0.2$ resuelva los siguientes ejercicios:
a)Encuentre las ecuaciones de movimiento en coordenadas cartesianas y polares. No es necesario resolver las ecuaciones en ambos sistemas, elija el que más le guste y con una transformación de coordenadas obtenga los resultados para el otro sistema de coordenado.
b)Utilice el método de Runge-Kutta de 4to orden para encontrar la solución a las ecuaciones de movimiento en el inciso anterior con condiciones iniciales $\vec{r}_0=(0,1)\ $ y $\ \vec{v}_0=(cos\theta,sen\theta)$, donde $\theta=\frac{n\pi}{6}$ con $n=0,1,...,11$ .
c)Genere una función que calcule le energía para distintas condiciones iniciales y detenga la simulación cuanto ésta sea menor al 1% de su valor inicial.
d)Dibuje las gráficas del momento angular para los casos anteriores del inciso b).
e)Tome de manera arbitraria una condición inicial y dibuje el espacio fase $(x,p_x)\ $ y $\ (\rho,p_\rho)$.
Solución¶
a)Por 2a Ley de Newton se tiene que $m\ddot{\vec{r}}=-k\vec{r}+\vec{f}_{fric}=\begin{cases}-k\vec{r} -\gamma\vec{v} &, \mbox{si } |\vec{v}|<1 \\-k\vec{r} -\mu\ v^{3/2} \hat{v} &, \mbox{si } |\vec{v}|>1 \end{cases}$. Entonces:
En Coordenadas Cartesianas:Tomamos que $\vec{r}(t)=(x(t),y(t))=x(t)\hat{i}+y(t)\hat{j}$,con $\hat{i}=(1,0), \hat{j}=(0,1)$. Con esto, la velocidad y la aceleración son símplemente: $\vec{v}=\dot{\vec{r}}=\dot{x}\ \hat{i}+\dot{y}\ \hat{j} \Rightarrow \vec{a}=\ddot{\vec{r}}=\ddot{x}\ \hat{i}+\ddot{y}\ \hat{j}$. Entonces, las ecuaciones de movimiento estan dadas por:
$$\ddot{x}_i=\begin{cases}-\frac{k}{m}x_i-\frac{\gamma}{m}\dot{x}_i &\text{, si }\ (\dot{x}_1^2+\dot{x}_2^2)^{1/2}<1 \\ -\frac{k}{m}x_i-\frac{\mu}{m}(\dot{x}_1^2+\dot{x}_2^2)^{1/4}\dot{x}_i &\text{, si }\ (\dot{x}_1^2+\dot{x}_2^2)^{1/2}>1 \end{cases}, \text{donde } i=1,2\ \text{ y }\ (x_1,x_2)=(x,y). $$Para resolver la ecuación hacemos el cambio de variable $\begin{cases}x_1=x(t)\\x_2=y(t) \\x_3=\dot{x}(t)\\x_4=\dot{y}(t) \end{cases} \Rightarrow \begin{cases}\dot{x_1}=x_3\\ \dot{x_2}=x_4 \\ \dot{x_3}=\ddot{x}(t)\\ \dot{x_4}=\ddot{y}(t) \end{cases}$.
Por tanto, el sistema se puede escribir como $\dot{\vec{x}}=\vec{g}(t,\vec{x})$, con $\vec{x}=(x_1,x_2,x_3,x_4)$ y $$\vec{g}=\begin{cases}(x_3,x_4,-\frac{k}{m}x_1-\frac{\gamma}{m}x_3,-\frac{k}{m}x_2-\frac{\gamma}{m}x_4) &\text{, si }\ (x_3^2+x_4^2)^{1/2}<1 \\(x_3,x_4,-\frac{k}{m}x_1-\frac{\mu}{m}(x_3^2+x_4^2)^{1/4}x_3,-\frac{k}{m}x_2-\frac{\mu}{m}(x_3^2+x_4^2)^{1/4}x_4) &\text{, si }\ (x_3^2+x_4^2)^{1/2}>1 \end{cases}$$.
En Coordendas Polares Planas:Tomamos que $\vec{r}(t)=r(t)\hat{r}(t)$,con $\hat{r}=(cos\theta(t),sen\theta(t))$. Con esto, la velocidad y la aceleración son: $\vec{v}=\dot{\vec{r}}=\dot{r}\ \hat{r} + r\dot{\theta}\ \hat{\theta} \Rightarrow \vec{a}=\ddot{\vec{r}}=(\ddot{r}-r\dot{\theta}^2)\ \hat{r}+(2\dot{r}\dot{\theta}+r\ddot{\theta})\ \hat{\theta}$, donde $\hat{\theta}=(-sen\theta(t),cos\theta(t))$. En este caso, tenemos una ecuación diferencial en cada dirección coordenada, pero las ecuaciones son distintas. Después de hacer todos los cálculos resulta que:
$$\begin{array}{rl}\text{para }\hat{r}:\quad & \begin{cases} r\dot{\theta}^2-\frac{k}{m}r -\frac{\gamma}{m}\dot{r} & \text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}<1 \\ r\dot{\theta}^2-\frac{k}{m}r -\frac{\mu}{m}(r^2\dot{\theta}^2++\dot{r}^2)^{1/4} \dot{r} &\text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}>1 \end{cases} \\ \text{para }\hat{\theta}:\quad & \begin{cases} -2\frac{\dot{r}\dot{\theta}}{r} -\frac{\gamma}{m}\dot{\theta} &\text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}<1 \\ -2\frac{\dot{r}\dot{\theta}}{r} -\frac{\mu}{m}(r^2\dot{\theta}^2+\dot{r}^2)^{1/4}\dot{\theta} & \text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}>1 \end{cases} \end{array} $$Para resolver la ecuación hacemos el cambio de variable $\begin{cases}r_1=r(t)\\r_2=\theta(t) \\r_3=\dot{r}(t)\\r_4=\dot{\theta}(t) \end{cases} \Rightarrow \begin{cases}\dot{r_1}=r_3\\ \dot{r_2}=r_4 \\ \dot{r_3}=\ddot{r}(t)\\ \dot{r_4}=\ddot{\theta}(t) \end{cases}$.
Por tanto, el sistema se puede escribir como $\dot{\vec{r}}=\vec{g}(t,\vec{r})$, con $\vec{r}=(r_1,r_2,r_3,r_4)$ y $$\vec{g}=\begin{cases}(r_3,r_4,r_1 r_4^2-\frac{k}{m}r_1-\frac{\gamma}{m}r_3,-2\frac{r_3 r_4}{r_1}-\frac{\gamma}{m}r_4) &\text{, si }\ (r_1^2 r_4^2+r_3^2)^{1/2}<1 \\(r_3,r_4,r_1 r_4^2 -\frac{k}{m}r_1-\frac{\mu}{m}(r_1^2 r_4^2+r_3^2)^{1/4}r_3,-2\frac{r_3 r_4}{r_1}-\frac{\mu}{m}(r_1^2 r_4^2+r_3^2)^{1/4}r_4) &\text{, si }\ (r_1^2 r_4^2+r_3^2)^{1/2}>1 \end{cases}$$.
b)Utilizamos el método de Runge-Kutta de 4o Orden. Para eso utilizaremos el código que se implementó en las tareas previas.
import numpy as np
import pylab as pl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
%matplotlib inline
def arg_rk4(funcion,x0,tlist,args=0,h=0.1):
"""
Funcion de Integracion por el Metodo de Runge-Kutta a 4o orden. Sus argumentos son: la funcion de la ecuacion diferencial,
el valor de la funcion al tiempo inicial x0,una lista con el tiempo inicial y el final tlist=[ti,tf] y
el paso de integracion como argumento opcional.
Puede tener una lista de parametros opcionales args que se pasan a la funcion de la ecuacion diferencial.
Regresa los tiempos y sus respectivos valores de la funcion como una tupla de arrays t,x.
Cabe notar que la funcion, el jacobiano y x0 pueden tomar valores vectoriales, los cuales se representan por arrays.
Tambien, la funcion debe tener ordenados los parámetros de la forma g(t,x).
"""
if (tlist[1]-tlist[0])<10.0*h: #cambia el valor de h para garantizar que se hagan almenos 10 iteraciones
h=(tlist[1]-tlist[0])*0.1
tiempos=np.arange(tlist[0],tlist[1],h) #regresa un array
x = np.zeros((len(tiempos),len(x0)))#regresa otro array del ancho de la dimension de x0 y el largo del array de tiempos)
x[0,:] = x0 #x[tiempo,coordenada],:==todas las coordenadas
if args == 0: #no hay argumentos opcionales dados
for i in xrange(0,len(tiempos)-1):
k1=funcion(tiempos[i],x[i,:])
k2=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k1)
k3=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k2)
k4=funcion(tiempos[i]+h,x[i,:]+h*k3)
x[i+1,:]=x[i,:]+(1/6.0)*h*(k1+2*k2+2*k3+k4)
else:
for i in xrange(0,len(tiempos)-1):
k1=funcion(tiempos[i],x[i,:],args)
k2=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k1,args)
k3=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k2,args)
k4=funcion(tiempos[i]+h,x[i,:]+h*k3,args)
x[i+1,:]=x[i,:]+(1/6.0)*h*(k1+2*k2+2*k3+k4)
return tiempos,x
#regresa una tupla, que la hace inmutable a diferencia de una lista y que, a diferencia de un array,
#cada entrada puede ser de distinto tipo (heterogéneas)
Ahora programamos las ecuaciones de movimiento en los distintos sistemas coordenados.
def mov_cartesianas(t,x_descartes,k=5.,m=1.,gamma=0.1,mu=0.2): #x_descartes=(x1,x2,x3,x4)=(x,y,vx,vy)
x=x_descartes #asumiremos que este es un array
v=np.sqrt(x[2]**2+x[3]**2)
if v<=1:
dx=[x[2],x[3],-k*x[0]/m-gamma*x[2]/m,-k*x[1]/m-gamma*x[3]/m]
else:
dx=[x[2],x[3],-k*x[0]/m-mu*np.sqrt(v)*x[2]/m,-k*x[1]/m-mu*np.sqrt(v)*x[3]/m]
return np.array(dx) #regresa un array
def mov_polares(t,x_polares,k=5.,m=1.,gamma=0.1,mu=0.2): #x_polares=(r1,r2,r3,r4)=(r,theta,vr,vtheta)
r=x_polares #asumiremos que este es un array
v=np.sqrt(r[0]**2*r[3]**2+r[2]**2)
if v<=1:
dr=[r[2],r[3],r[0]*r[3]**2-k*r[0]/m-gamma*r[2]/m,-2*r[2]*r[3]/r[0]-gamma*r[3]/m]
else:
dr=[r[2],r[3],r[0]*r[3]**2-k*r[0]/m-mu*np.sqrt(v)*r[2]/m,-2*r[2]*r[3]/r[0]-mu*np.sqrt(v)*r[3]/m]
return np.array(dr) #regresa un array
Ahora, utilizamos el metodo de Runge-Kutta para integrar las ecuaciones de movimiento, utilizando condiciones iniciales $\vec{r}_0=(0,1)\ $ y $\ \vec{v}_0=(cos\theta,sen\theta)$, donde $\theta=\frac{n\pi}{6}$ con $n=0,1,...,11$ .Además de simplemente dibujar todas las soluciones en una o varias gráficas optaremos por animar la solución a fin de que el resultado sea más ilustrativo.
Coordenadas Cartesianas
paso=0.01
tf=10
x0,y0=0,1 #px=(x,y,vx,vy)
fig = plt.figure(figsize=(5,5))
colormap = plt.cm.gist_rainbow
plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, 12)])
sols_cart=[]
tiempos_cart=np.arange(0,tf,paso)
for n in xrange(12):
theta=np.pi*n/6
vx0=np.cos(theta)
vy0=np.sin(theta)
px=[x0,y0,vx0,vy0]
tx,x = arg_rk4(mov_cartesianas, px, [0,tf],h=paso)
sols_cart.append(x)
plt.plot(x[:,0],x[:,1],lw="2")
plt.plot(0,1,marker='o',color="black",markersize=10)
plt.title('Solucion en Coordendas Cartesianas')
plt.show()
from JSAnimation import IPython_display
from matplotlib import animation
# Se define el ambiente en el que queremos hacer la animación
fig = plt.figure()
ax = plt.axes(xlim=(-2, 2), ylim=(-2, 2))
lines = [plt.plot([], [],lw=2)[0] for j in range(12)]
points = [plt.plot([], [],marker='o',color="black")[0] for j in range(12)]
T = 24*tf # periodo para generar 2*T cuadros
# Funcion para inicializar cada cuadro de la animacion
def init():
for line in lines:
line.set_data([], [])
for point in points:
point.set_data([],[])
return lines,points
# Esta funcion se llama de manera secuencial para cada elemento k.
step=len(tiempos_cart)/(2*T)
def animate(k):
for j,line in enumerate(lines):
xdata=sols_cart[j][:,0]
ydata=sols_cart[j][:,1]
line.set_data(xdata[1:(step*k)],ydata[1:(step*k)])#grafica una linea con los puntos previos
for j,point in enumerate(points):
xdata=sols_cart[j][:,0]
ydata=sols_cart[j][:,1]
point.set_data([xdata[step*k]],ydata[step*k])#grafica sólo el punto actual
return lines,points
# Se llama a la animacion. blit=True es para que solo se dibije las partes de la imagen que tienen cambios.
animation.FuncAnimation(fig, animate, init_func=init,frames=2*T, interval=20, blit=True)